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ABSTRACT 

We present a method for computing the 6-dimensional coarse-grained phase-space den- 
sity /(x, v) in an N-body system, and derive its distribution function v{f). The method is 
based on Delaunay tessellation, where v{f) is obtained with an effective fixed smoothing 
window over a wide / range. The errors are estimated, and v{f) is found to be insensitive 
to the sampling resolution or the simulation technique. We find that in gravitationally re- 
laxed haloes built by hierarchical clustering, v{f) is well approximated by a robust power 
law, v{f) oc j-2.5±o.05^ Qygj. jjjore than 4 decades in /, from its virial level to the numerical 
resolution limit. This is tested to be valid in the ACDM cosmology for haloes with masses 
10^ — indicating insensitivity to the slope of the initial fluctuation power spectrum. 

By mapping the phase-space density in position space, we find that the high-/ end of v{f ) 
is dominated by the "cold" subhaloes rather than the parent-halo central region and its global 
spherical profile. The value of / in subhaloes near the virial radius is typically > 100 times 
higher than its value at the halo centre, and it decreases gradually from outside in toward its 
value at the halo centre. This seems to reflect phase mixing due to mergers and tidal effects 
involving puffing up and heating. The phase-space density can thus provide a sensitive tool 
for studying the evolution of subhaloes during the hierarchical buildup of haloes. It remains 
to be understood why the evolved substructure adds up to the actual universal power law of 
v{f) cx f~^'^. It seems that this behaviour results from the hierarchical clustering process 
and is not a general result of violent relaxation. 

Key words: cosmology: theory — dark matter — galaxies: dwarfs — galaxies: formation — 
galaxies: haloes — galaxies: interactions — galaxies: kinematics and dynamics 



1 INTRODUCTION 

The standard paradigm assumes that dark-matter haloes are the ba- 
sic entities in which luminous galaxies form and live. The haloes 
dominate the gravitational potential over a wide range of radii and 
they have a crucial role in determining the galaxy properties. While 
many of the systematic features of halo structure and kinematics 
have been revealed by A^-body simulations, the origin of these fea- 
tures is still an open theoretical puzzle, despite the fact that they are 
due to simple Newtonian gravity. 

The halo density profile p(r) is an example for such a 
puzzle. It is found in the simulations to have a robust non- 
power-law shape (which we refer to in general as "NFW"), 
with a log slope —3 near the virial radius, flattening gradu- 
ally toward —1 at about 1% of the v irial radius, and perhaps 
flattening further at smal l er radii, (e.g . |Navano,Jrenk^WhiteJ 
| l997t lMoore etall|l998t iKlvpin et al 1 1200 it IPower et al 1 120031 



file is insensitive (or at most weakly sensitive) to parameters of 
the cosmological model and the initial fluctuation power spec- 
Cole c& Lacev 1996; Navarro, Frenk& White 1997 
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200C: iRicotti l l200i 


Colin et al 


l200i .Navarro et al.i2004) indicating that its origin 


IS due to a 



robust relaxation process rather than specific initia l conditions. 
In particular, violent relaxation iLvnden-B ell 1967) may be in- 
volved in shaping up the density profile. In addition, secondary 
infall may be important in the outer regions Eun n and Gqty 

ichl 



'Dekel, Kowitt & Shaham' 
Hoffman c& Shaham 1985; 
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1981; Fillmore and Gol dreich I 
White & Zaritskv 1992J while 



Hav^liie^m2004 , and references therein). This density pro- 



some argue that resonances may have a r ole in the irmer region s 
( Weinberg &Katz 2002; but see iValenzuela & Klvpinl l2003h . 
Nevertheless, there is no clear understanding for why the haloes 
actually pick up the particular density profile they have. 

The properties of the velocity dispersion tensor is another the- 
oretical puzzle. The haloes tend to be rather round, with a velocity 
dispersion profile that is slightly rising at small radii and slightly 
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falling at large radii but is rather flat overall taussetal 11999 Jd) . 
The profile of the anisotropy parameter /3(r), which is a measure of 
radial versus tangential velocities, indicates near isotropy at small 
radii, which gra dually dev elops into more radial orbits at large radii 
jCole&LacevU 996: H uss et alll999JbllColm et all2000j) . For 
a spherical system in equilibrium, the a(r) and /3(r) are related to 
p(r) via the Jeans equation, but it is not at all clear why the haloes 
in the simulations choose the characteristic (j(r) or /3(r) that they 
actually do. 

An interesting at tempt to address the origi n of the halo pro- 
file has been made by iTavlor & Navarro I i200lh . who measured a 
poor-man phase-space density profile by /tn('') = p{r)/a{r)'^, 
and found that it displays an approximate power-law behaviour, 
/tn oc r~^ *^, over more than two decades in r. Using the spheri- 
symmetric Jeans equation, they showed that this power law permits 
a whole family of density profiles, whose limiting case is a profile 
similar to NFW, which asymptotically approaches a slope —0.75 
as r — > 0. The general power-law shape of /tn (r) is confirmed in 
the simulated haloes described below. This scale-free behaviour of 
/tn (r) is intriguing, and it motivates further studies of halo struc- 
ture by means of phase-space density. 

Simulations of the ACDM cosmology also reveal a roughly 
self-similar hierarchical clustering process, where smaller build- 
ing blocks accrete and merge into bigger haloes. At every moment, 
every halo contains a substructure of subhaloes on top of a smooth 
halo component that has been tidally stripped from an earlier gener- 
ation of substructure. Some of the important dynamical processes 
involved in this hierarchical halo buildup are understood qualita- 
tively. These include, for example, the dynamical friction which 
governs the decay of the satellite orbits, the tidal stripping of sub- 
haloes due to the host halo potential, and the mergers and flyby 
interactions of the subhaloes among themselves. However, a com- 
plete theoretical understanding of how these processes work in de- 
tail, and how they combine to produce the halo structure and kine- 
matics, is lacking. 

Attempts have been made to explain an inner density 
cusp using toy models of dynamical stripping and tidal ef- 
fects during the halo buildup by mergers jSver & White 
1998'; 'Nusser& ShethI Il999t iDekel. Devor & Hetzroni I 120031: 
Dekel et al. 2003). However, a similar halo density profile seems 
to be produced also in some of th e simulations where s ubstructure 
has be en artificially suppressed jMoo^^^J^99bljHus^^lJ 
| l999j:lAv ila-Reese et al. 200 l|jBullock. Kravtsov & Colin l200a 
lAlya rez^hapiro & Martel 200?), indicating that the process re- 
sponsible for the origin of this density profile might be a somewhat 
more general feature of gravity and not unique to the merger sce- 
nario. 

The issue of halo substructure has become timely both be- 
cause of its relevance to observations and its implications on other 
major issues in galaxy formation. Tidal tails and streams associ- 
ated with dwarf satellite gala xies have been observed in the haloes 
of the Milky Way and M31 ^I bata et al |[T99l l200lJbl) . and they 
start to allow detailed modelling of the halo history through the 
satellite orbits. Gravitational-lensing observations provide prelim- 
inary indications for the presence of substructure in haloes at the 
high level predicted by t he dissipationless ACDM scenario (e.g. 
iDalal & Ko chanek '2002^ In contrast, the observed number den- 
sity of dwarf galaxies seems to be significantly low er, thus pos- 
ing a "m issing dwarf problem" Klvcin et al. ( 1999b"): lMoore et all 
h999iJ) . Also, the " angular-mo mentum problem" of disk galax- 
ies (e.g. iNavarro & Steinmetz] l200d iBullock et aD I^OOlb) is 
probably associated with the evolution of substructure in haloes 



jMaller & Dekel l2002tlM'aller. Dekel & Somerville l2002h . While 
the dwarf and angular-momentum problems necessarily involve 
baryonic processes, understanding the gravitational evolution of 
substructure is clearly a key for solving them. 

In order to better understand the origin the various aspects 
of halo structure and buildup mentioned above, we make here a 
first attempt at addressing directly and in some detail the phase- 
space structure of dark-matter haloes. The fundamental quantity 
in the dynamical evolution of gravitating systems is the full, six- 
dimensional, coarse-grained, phase-space density /(x,v), which 
intimately relates to the underlying Vlasov equation and lies be- 
hind any relaxation proces s that may give rise to the virialized halo 
structure iBinnev & Tr emaine 1987, chapter 4). 

Ideally, one would have liked to compute it free of assump- 
tions concerning spherical symmetry, isotropy, or any kind of equi- 
librium. However, computing densities in a six-dimensional space 
is a non-trivial challenge which requires simulations of a very 
broad dynamical range. The state-of-the-art N-body simulations, 
with more than million particles per halo, allow an attempt of this 
sort for the first time. We describe below a successful algorithm for 
measuring /(x, v), and study its relevant properties including the 
associated systematic and random uncertainties. We then apply this 
algorithm to simulated virialized haloes in the ACDM cosmology. 

We report in this paper two surprising new results. First, we 
discover that the volume distribution function of the phase-space 
density, «(/), displays a universal scale-invariant power-law shape, 
valid in all virialized haloes that form by hierarchical clustering. 
Second, we realise that this power law is not directly related to the 
overall density profile, but is rather driven by the halo substruc- 
ture. This implies that the phase-space density provides a useful 
tool for studying the hierarchical buildup of dark-matter haloes and 
the evolution of substructure in them. 

In Sj2|we introduce /(x,v) and v{f). In S|3|we describe the 
method of computing / and «(/) from an A^-body halo, and sum- 
marise its properties and the associated errors, which are addressed 
in more detail in Appendix]^ In f|4|we present the universal power- 
law shape of v{f) based on several different simulations, and 
demonstrate its robustness to the mass scale and simulation tech- 
nique. In S|5|we display maps of phase-space density and show that 
the high-/ contributions to v{f) come from substructures within 
the parent halo. In f|6|we summarise and discuss our results and 
future work. 



2 THE PDF OF PHASE-SPACE DENSITY: «(/) 

This is a more technical and elaborate introductory section, aimed 
at introducing the concepts and nomenclature relevant for the anal- 
ysis of this paper. 

2.1 Definitions 

The state of a collisionless system is completely determined by the 
fine-grained phase-space density function /(x, v,t), which mea- 
sures the mass contained in an infinitesimal phase-space patch of 
volume dxdv, located at (x, v). The evolution of /(x, v, t) is gov- 
erned by the Vlasov equation, 

at/ + v- Vx/- Vx$- Vv/ = o, (I) 

with "I>(x) the gravitational potential, related self-consistently to 
/(x, v) by the Poisson's integral 
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$(x) 



G dx'dv 



/(x',v) 



(2) 



It is therefore sensible to assume that a true understanding 
of the nature of self-gravitating collisionless systems must involve 
/(x, v) as a primary ingredient. However, being a function of six 
variables, /(x, v) is hard to deal with computationally. It turns out 
that there is a much simpler function, v{f), which is intimately re- 
lated to /(x, v), but is much easier to handle both analytically and 
numerically. 

Before defining v{f), it is worth recalling that the Vlasov 
equation preserves the phase-space density along particle trajec- 
tories: It implies eq. Q implies that 



d_ 

di 



/[x(t),v(f)] = 0, 



(3) 



where [x(t), v(t)] is the trajectory of a phase-space element. This 
conservation of phase-space density along trajectories can also be 
viewed as a conservation of phase-space volume. In order to il- 
lustrate what this means, consider any smooth function C{x), and 
define the integral 



I{t) = JdxdvC[f{x,v,t)] . 
Using the Vlasov equation we find that 



dxdvC'(/)9t/ 



J dxdv C'if) [v ■ Vx/ - Vx<& • Vv/] 



(4) 

(5) 
(6) 



= -J dxdv [v ■ VxC(/) - Vx$ ■ VvC(/)] = . (7) 

The last equality follows by integrating the first term over x and the 
second term over v, and assuming that / ^ as x, v — » oo. This 
implies that / is conserved under the dynamics. 

Now, if C(/) is the Dirac delta function^ S{f — /o), then the 
integral / becomes an /-dependent function, which we refer to as 
the PDF v{f ): 



v{f = fo) 



Jdxdv5[f{yL,v,t)- fo] 



(8) 



From this definition it is clear that v{f)df is the volume of phase- 
space occupied by phase-space elements whose density lies in the 
range (/, / + df). The conservation of v{f) implied by Eq. 
means that this volume is conserved under the dynamics. The fine- 
grained v{f) [i.e., the v{f) which is calculated using the fine- 
grained /(x, v)] should therefore be viewed as a signature of the 
system, which remains the same throughout its evolution. 



2.2 Coarse-Grained Density and Mixing 

One might have thought that the conservation of v{f) poses severe 
constrains on the evolution of /(x, v,t), but this is not the case 
due to the effects of mixing. In the course of evolution, phase-space 
patches of high / are stretched and spiral into regions with low /. 
As the stretching continues, the patches become thiner and thiner, 
and as a result /(x, v) varies over increasingly smaller scales. Very 
soon, one can no longer measure /(x, v), but instead measure an 



^ Strictly speaking, the Dirac delta function is not a smooth function, but it 
can be approximated by a series of smooth functions, each obeying Eq. 171 



average of it over some finite volume. This average is often re- 
ferred to as the "coarse-grained" phase-space density, as opposed to 
the original "fine-grained" density /(x, v). We denote the coarse- 
grained quantities by / and v{f), but in subsequent sections we 
may omit the bar and simply refer to them as / and v{f). The im- 
portant point to realise is that the coarse-grained v{f) is no longer 
conserved. For example, given initially a volume V filled with 
phase-space density /a and a similar volume filled with density 
/s, one ends after mixing with a volume 2V filled with an average 
density (/a + /s)/2. 

When a large fraction of the mass is added to halo, e.g. by 
collapse or a major merger, rapid global fluctuations of the gravita- 
tional potential re-distribute the energies of each phase-space ele- 
ment, and lead to very strong mixing across large scales, resulting 
in variations in /. After a few global dynamical times, the potential 
fluctuations fade away and / stabilises. The fine-grained /, on the 
other hand, continues to spiral and stretch indefinitely, on smaller 
and smaller scales, which no longer affect /. After a while, the / 
can be viewed as the "physical" phase-space density of the system, 
since the microscopic fluctuations of / can no longer be measured, 
and, in particular, no longer affect the gravitational potential. It is 
therefore the coarse-grained / which, according to Jeans' theorem, 
can be written as a function of he invariants, the energy E and the 
angular momen tum L. This relaxat ion process is referred to as vio- 
lent relaxation i Lvnden-Bellll967h . 

What is the coarse-grained v( f) of a relaxed system, and how 
is it related to the constant fine-grained v[f) of that system? One 
may identify the fine-grained v{f) with the coarse-grained v{f) at 
an early time when the initial / is smooth enough not to vary over 
scales smaller than the averaging scale associated with /. Then, the 
question can be rephrased in terms of the relation between the final 
v(f) and the early v{f). One obvious constraint on the final /, as 
an average of /, is that their maximum values should obey /max ^ 
/max. In particular, if the early v{f) [namely v(/)] vanishes for 
/ > /max, then so does t he final t)(/). 

The mixing theorem frremaine. Henon & Lvnden -Bell 198a) 
specifies several additional constraints on the final v{f), which 
arise from the fact that it has originated from a given «(/). The 
cumulative distribution function is defined by 



dxve[/(x,v)-/o] 



v{f)df , 



(9) 



fa 



where 6[-) is the Heaviside step function. As a monotonic func- 
tion, V{f) can be inverted to yield f{V), the reduced distribution 
function. The associated cumulative mass function M{f), which 
measures the mass in phase-space regions where /(x, v) < fo, is 
then defined by 

M(/o) = y"dxdv/(x,v)e[/(x,v)-/o] 

f'vindf. 

fa 

Substituting / — f{V) inside M{f), one obtains the function 
AI{V), which can also be written as 



(10) 
(11) 



M{V) = / /(V') dV' 



(12) 



If Af(V) and M'{V) denote the functions arising from the density 
functions /(x, v) and /'(x, v) respectively, then the mixing theo- 
rem states that the final (coarse-grained) /' is an average of the ini- 
tial (coarse- or fine-grained) / if and only if M'(V) ^ M{V) for 
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every V . This, in turn, implies an implicit constraint on the evolu- 
tion of the coarse-grained v{f). This is one of the very few rigorous 
results concerning the evolution of v{f), but its strength is rather 
limited because it provides only an integro-differential inequality 
that constrains «(/). 

The mixing theorem c an be also stated directly in terms of 
v{f), as shown by Mathur Hl988l) . In his formulation, a necessary 
condition for the final (coarse-grained) v{f) evolved from the ini- 
tial v{f) is that^ 

v{f)=<f) + ^P{f), (13) 

with P{f) being some continuous function (not to be confused 
with a probability distribution) such that P{f) ^ for all / and 
P(/) ^ for / ±oo. 

In the CDM scenario, the dark matter is initially cold - with 
vanishing velocity dispersion. This means a one-to-one correspon- 
dence between the particle positions and their velocities, v = 
v(x), and the initial (fine-grained) phase-space density can be writ- 
ten as 



/o(x,v) =p(x)5^[v-v(x)] 



(14) 



Thus, in the beginning, the system can be described by the evolu- 
tion of its density and velocity fields, often approximated by the 
Zel'dovich approximation. This singular nature of the initial phase- 
space density, together with its conservation, plays a crucial role in 
the formation of structures, and in particular in the abundance of 
dwar f haloes, and in the formation of cusps in dark-matter haloes 
(e.g. IColm et a r2000b"). 

Physically, the phase-space density is never really infinite, as 
it is defined by a finite number of dark-matter particles per given 
phase-space volume. Numerically, it is the finite mass-resolution of 
the A'^-body simulation that prohibits the phase-space density from 
being infinite. Therefore, the initial phase-space density is expected 
to be extremely high in regions where v ~ v(x) and vanishingly 
small elsewhere. Accordingly, v{f ) will have contributions from 
a narrow range of extremely high values of / and from a narrow 
range near / = 0. As the system evolves and mixing takes place, 
phase-space densities from both regions mix and give rise to in- 
termediate values of /, widening the distribution v{f ) under the 
constraints of the mixing theorem. 



2.3 A Relation Between v{f) and p(r) 

In general, different systems may have the same v{f). However, 
if the system is spherically symmetric and stationary, such that / 
is a function of the energy alone, /(x, v) = /(e), then there is a 
unique relation between v{f), /(e) and p{r). We shall see in Sj4] 
that the v{f) measured in A'^-body haloes is actually different from 
the v{f) one would have predicted from the halo p(r) using the 
relation assuming /(e). 

For a spherical system the relation between p(r) and v{f) can 
be derived as follows. Assume that /(x,v) = /[e(x, v)], where 
e(x, v) = /2 + "I>(r), and define the density-of-states function 



5(e) 



The quantity g{e)de measures how much phase-space volume is 
occupied by phase-space elements with energy in the interval 
(e, e + de). When the system is spherical, g{e) can be written as 
i Binney &Tremaine 19 87., Eq. 4- 157b), 

gie) = (ATrf [ ^ \^ J 2[e - <!>is)] ds , (16) 
Jo 

with r-(e) the inverse function of the gravitational potential. The 
overall phase-space volume occupied by energy levels below e is 
thus given by 



V 



g{e')de 



:(* 



n)' r'".-^{2[e-$(s)]}'^' ds. 
Jo 



(17) 



(18) 



Therefore, if V{f) is the cumulative v{f) defined in Eq. j9j, and 
/(x, v) is a function of the energy alone, we obtain an integro- 
differential equation for /(e): 



■(<=) 



(19) 



This has to be supplemented by the equation that connects c[>(r) to 
/(e) jBinnev & Tremainell987i eq. 4-104), 



dr J 



47rGp(r) 



(20) 



= {47:fG /(e)v/2($-e)de. (21) 

J*(r) 

In general, these equations can be solved numerically, e.g., using a 
Picard iteration scheme. 

In the asymptotic regime, the above equations can be solved 
analytically. If p(r) oc r~°' as r — + 0, then, for e $(0), one can 
show that /(e) and g{e) have the scale-invariant forms 



(22) 
(23) 

(24) 



dxdv 5[e(x, v) 



(15) 



/(e) cx [e-0(O)] 2(2-) , 
g{e) oc [e^<l>{0)]^^ . 

Note that the derivative of Eq. <19> with respect to e yields 

rr/ M ^(e) 

v[.m] = ^ . 

By plugging the scale-invariant solutions of eq. <23> into eq. <24t . 
one finds that for / oo, the PDF is also a power law, v{f) oc 
/"'', with 

18 -4q 18-6/3 
P = or a = — ; • (25) 

We learn that the density slopes in the range a ^ 2, relevant 
for the inner regions of haloes where / is high, correspond to a 
narrow range of (3 values, 3^/3^ 2.5. In this case of power- 
law profiles, l3 = 2.5 corresponds to the singular isothermal sphere 
a — 2, P — 2.8 corresponds to an a = 1 cusp, and /? = 3 
corresponds to a flat core, q = 0. 

Eq. 03 can also be obtained from simpler dimensional argu- 
ments using the virial theorem. 



This is the original statement in However, we believe 

that by adding the (trivial) condition that v{f) ^ for every /, it becomes 
also a sufficient condition. 



3 MEASURING vif) IN AN iV-BODY SYSTEM 

We wish to measure the phase-space density /(x, v) of a system 
represented by A'^ particles of mass m each. One straightforward 
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approach would have been to divide phase-space into a large num- 
ber of Cartesian cells of equal volume V each. If cell j contains 
Nj particles, then the average phase-space density in it can be esti- 
mated as fj — Njin/V. However, this approach is impractical in 
a six-dimensional space. A moderate resolution of 100 bins along 
each axis would require lO'^^ cells, which is beyond the capacity of 
any present-day computer. Moreover, for the computed / to have 
any statistical meaning, there should be at least 100 particles in 
each cell (for a relative Poisson error of ~ 10%), which requires 
a total of 10^* particles - about 5 orders of magnitude beyond the 
number of particles in today's largest simulations. 

A possible way to overcome this problem is by using an adap- 
tive grid, where the cells vary in size to properly allow high res- 
olution in high-density regions and low resolution in low-density 
regions. One can simply create an adaptive grid by dividing each 
Cartesian cell which contains more particles than some prescribed 
threshold TV* into two or more sub-cells. This division can be done 
recursively until all cells contain N* particles or less. 

Even more effective would be to vary the shape of the cells 
as well, allowing them to adapt more efficiently to the geometry 
of the underlying distribution and thus increase the effective res- 
olution. A particularly robust method of this type is the Delau- 
nay Tessellation Field Estimator (DTFE), which has already been 
impleme nted in three dimensions in a cosmological context (e.g. 
[b ernarde au & van de Wevgaertll99d : ISchaap & van de Wevsaert I 
I2OOO . and references therein). We use this method here to measure 
/ in six dimensions. 



3.1 Constructing a Delaunay Tessellation 

A tessellation is the division of i?'' space into a complete cover- 
ing of mutually disjoi nt convex polygons. The Delaunay tessella- 
tion (Delaunay 1 11 934) is defined for a sample of A'^ points as fol- 
lows. The Delaunay cells that construct the tessellation are the d- 
dimensional polyhedrons made by connecting every set of d + 1 
points whose circumsphere [the (d — 1) -dimensional sphere that 
passes through all of them] does not encompass any other point 
from the sample. One immediate advantage of the Delaunay tessel- 
lation is that it is parameter-free, and it completely adapts itself to 
the un derlying distribution of points. Sph aap & van de Wevgaert 
i2OO0l) have demonstrated the superiority of the DTFE over more 
conventional methods for estimating the density in 3D (methods 
like cloud-in-cell or the smoothing kernel used in SPH simula- 
tions). One may assume that the same holds for the 6D case. 

In 2D, the Delaunay cells are triangles, as illustrated in Fig.Q 
In 3D, the Delaunay cells are tetrahedrons. In the six-dimensional 
phase space, the Delaunay cells are six-dimensional polyhedrons, 
each defined by 7 particles. 

Figure |2| shows the Delaunay tessellation of an uneven distri- 
bution of points in the 2D plane. It demonstrates the obvious adap- 
tive nature of the tessellation: regions with high density of particles 
are covered by small triangles, whereas regions with low density 
are covered with larger triangles. 

Constructing the Delaunay tessellation of A'^ lO" par- 
ticles in a six-dimensional space is not a straightforward task. 
Out of the many algorithms that exist in the literature (e.g., 
Su & Drvsdale 1995, and references therein), we followed 
Berna rdeau & van de Wevgaert ( 19 96) and [van de Wevgaert 
^^94!) in picking the algorithm by iTanemura. Oeawa & Oeital 
il983h . We included some of the programming adjustment by van 
de Weygaert, e.g., the use of fc-trees to speed up searching, and 
converted the code from three dimensions to six dimensions. 




Figure 1. The Delaunay tessellation of 8 points in two dimensions. Shown 
for example is a circumcircle of one Delaunay cell, which, by construction, 
does not encompass any other point. 




Figure 2. An illustration of the adaptive nature of the Delaunay tessellation. 
Left: an uneven distribution of points in 2D. Right: the resulting Delaunay 
tessellation. 



A full account of the algorithm, the numerical implementa- 
tion and code performance, is provided elsewhere (Arad 2004, in 
preparation). Here we just mention briefly that for a halo of ~ 10® 
particles the C code runs for about a week on a common PC with 
CPU ~2GHz and internal RAM ~1GB. The resulting tessellation 
is made of ~ 10^ Delaunay cells, where a typical particle is sur- 
rounded by ~ 7, 000 cells involving ~ 200 neighbouring particles. 



3.2 Recovering / and v{f) from tlie Tessellation 

Once the Delaunay tessellation is constructed in phase space, we 
use it to estimate /(x, v), g eneralising the method implemented in 
3D bv lSchaap & van de Weveaertl J200CI) . First we estimate fi for 
each particle i. We define a macro cell by joining all its surrounding 
Delaunay cells: 



(26) 



where {Dj} is the set of all Delaunay cells and {D,,. } is the subset 
of cells which contain the particle i as one of their vertices. Figure|3| 
shows such a macro cell in 2D. The estimated phase-space density 
at point i should be inversely proportional to the volume of the 
macro cell, 1 /| Wi \ . The proportionality factor must be greater than 
unity because the Wi cells of different particles partly overlap, and 
one can show that it should in fact be d + 1 in order to preserve the 
total mass of the system. In phase-space we therefore define 



(27) 



In order to estimate / at any general point (x, y), one might 
have used linear interpolation based on the 7 vertices of the cor- 
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Figure 3. An illustration of a macro cell Wt (grey area), the union of the 
Delaunay cells that surround the particle i. 



responding Delaunay cell D^. This means expressing /(x,v) as 
a linear function of the 6D vector uo = (x, v), namely fiuo) = 
A,j -1^ + B^, with Ai, and Bi, a 6D vector and a scalar whose seven 
unknown values are found by equating /(x, v) with fi at the seven 
vertices of the cell. However, with 10*^ particles involving ~ 10® 
Delaunay cells the linear interpolation is CPU-expensive and prac- 
tically impossible. Instead, we perform a zero-order interpolation 
where /(x, v) is constant inside each Delaunay cell - the average 
of the fi values at the seven vertices of D^, : 

1 



This has the advantage that any integral of the form 



/ = y^xdv * [/(x, v)] 



(28) 



(29) 



with ^'(■) some arbitrary function, can be easily estimated by the 



(30) 



In particular, the desired v{f) can be found by first computing 
its cumulative counterpart V(f), 



V{fo) = / vU')df= / dxdv 

J fo J/(x,v)>/o 

and then taking its derivative, 

dV{f) 



(31) 
(32) 



df 



In what follows, we often denote the DTFE-measured / and 
v(f) by /del and Udci in order to distinguish them from the exact 
quantities. 

3.3 Error Estimate 

When we measure the v{f ) of a cosmological A'^-body system us- 
ing the Delaunay tessellation, we should expect two types of errors. 
First are the errors in the underlying / associated with errors in the 
numerical simulation itself, such as errors due to two-body relax- 
ation effects, force estimation, time integration and so on. A way to 
estimate these errors is by re-simulating the same system with dif- 
ferent codes and with different sets of numerical parameters gov- 
erning the mass resolution and the force resolution. A systematic 
testing of this sort will be reported in an associated paper (Arad, 



Dekel & Stoehr, in preparation). In the current paper, we make a 
preliminary comparison of several different haloes simulated with 
different resolutions and with different codes. We find that all the 
virialized haloes tested recover almost the same shape of v{f). 

The second type of errors originate from the fact that we try 
to estimate a smooth /(x, v) from a finite set of particles using 
the Delaunay tessellation. Here we may encounter both statistical 
and systematic errors. Some of these errors would decrease as the 
number of particles is increased, whereas some are an inherent part 
of the DTFE independent of the mass resolution. In order to ob- 
tain a simple understanding of the nature of the DTFE errors, we 
first use a simple statistical model based on a Voronoi tessellation, 
which resembles the DTFE but lends itself more easily to analytical 
treatment. We then evaluate the DTFE errors empirically by mea- 
suring /dci(x,v) in synthetic systems were the particles represent 
a known /(x, v), and compare the results to the error-model pre- 
dictions. Since this error analysis is somewhat detached from the 
main focus of this paper, we describe it in detail in Appendix IXI 
Our wisdom regarding the DTFE properties and uncertainties can 
be summarised in three points as follows: 

(i) The measured /dd at each particle is chosen, to a good 
extent, from a probability distribution of the random variable 
/del //true. lu a typical realisation with 10** particles, the width of 
this distribution is about one decade, defining the range of fluctua- 
tion of /del about /true. As TV ^ oo, the shape of the distribution 
approaches an asymptotic limit p(/dei//true) with & finite width 
of about one quarter of a decade, as deduced from the DTFE of a 
Poisson distribution. Therefore, also in the infinite limit, the DTFE 
is expected to produce local fluctuations. 

(ii) The measured Vdei[f) can be viewed as a convolution of 
the true «(/) and a. fixed window fiinction, which is related to the 
probability distribution p(-) [Eq. <A5H : 



""dclif = fo) 



if).fo Pifo/f) df . 



(34) 



For distributions where «(/) is close to a power-law, the difference 
between t;dei(/) and i'true(/) is negligible over a large range of 
scales, as demonstrated in Figs. lA2l and IA3I This is a very useful 
feature of the DTFE-measured v{f). 

(iii) The relative statistical error in Wdei (/) is proportional to 
and can be approximated by [Eq. JA18H 



A(/)=c 



/(Vdel(/)) 



1/2 



(35) 



(33) with Vdei(/) = JJ°Vdei{f') df and m = M/N. The constant c 
is of order unity, and can be calibrated by comparing Eq. <35> to 
the actual error in lower-resolution measurements. In practise, this 
means that when TV ^ lO" or more, the statistical error in Udei(/) 
is negligible for a very wide range of /. Moreover, in regions where 
there are large statistical errors, they are likely to be overwhelmed 
by systematic errors. 



4 A UNIVERSAL SCALE-FREE v{f) 

We have analysed the v{f) of several different haloes, in three dif- 
ferent mass ranges, simulated within the ACDM cosmology with 
two different TV-body codes. 

In order to calculate the v{f) of a given halo, we find the halo 
centre using a simple max-density algorithm, and extract all parti- 
cles which lie within a distance R from its centre. The max-density 
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algorithm is based on an iterative counting in cells: at each itera- 
tion, space is divided into 8 equal cubic cells and the densest cell is 
chosen for the following iteration. The iterative process stops when 
the densest cell contains no more than 500 particles, and the cen- 
tre of that cell is defined as the centre of the halo. The radius R is 
typically set to be ~ 10 — 20% larger than the virial radius. We 
compromise on analysing ~ 10^ particles, such that we explore 
a significant dynamical range while the computation can be com- 
pleted in a few days. 

To estimate the typical phase-space density at the outer regions 
of the halo, we define f^ir by 



fvi 



Pvi 



t3/2 , 



(36) 



which is the average space-space density that one would measure 
for a halo of constant density p^ir and a Maxwellian velocity distri- 
bution with a velocity dispersion avir- For a given halo in a given 
cosmology, the mean virial quantities p^ir and a^ir are defined 
using the virial theorem and the top-hat model relevant to that cos- 
mology. 

As one crude estimate of the upper limit on / for which the 
measurement of v{f) is reliable, we evaluate in each halo the / 
value, f2o%, below which the statistical error in v{f) due to the 
DTFE procedure is below 20%. For that we use Eq. 03, cah- 
brated by measurements with 10^ particles. This statistical error is 
expected to be practically negligible in the range /vir < f < .f20% ■ 

4.1 A^-body Simulations on Different Scales 

The results described in this paper are based on three different cos- 
mological simulatio ns. Two of them used th e Adaptive Refinement 
Tree (ART) code iKravtsov. Klvpin & Khokhloy .,.1997.). and the 
third used the Tree Pa rticle Mesh (TPM) code jSode et alll2000l: 
iBode & Ostriker 1200 3*). In all the simulations, the assumed cosmo- 
logical model is the standard ACDM with Sim = 0.3, f^A — 0.7 
and h = 0.7 today. 

The ART simulations were done in periodic boxes of sizes 
L = 1 ft~^Mpc and L = 25 h'^Mpc, whereas the TPM simu- 
lation was done in a box of L = 320 h^^Mpc. We denote these 
three simulations by LI, L25 and L320 respectively. Three haloes 
were analysed from each of these simulations, with masses corre- 
sponding to to dwarf galaxies (10® — 10^° Mq), normal galaxies 
(~ 10^^ M0) and clusters (~ W^'' Mq) respectively. Global prop- 
erties of of the simulations and the haloes are given in table l4m 
The v{f) curves for these haloe s are shown i n Fig. ID 

The LI simulation is from lColm et all i2003l) . It is analysed 
at z = 2.33 in physical coordinates. The three haloes studied, 
labelled A, B, C, are the largest haloes in the snapshot. Their 
virial radii refer to a mean overdensity of A = 183, as appro- 
priate fb£_the_given cosmology and redshift. The L25 simulation 
is bv lKlvpinetal] <200ll) . The haloes are denoted B, C, D fol- 
lowing the notation in the simulation paper. The L320 simulation is 
by|Wambs2anss et al (2004); Weller et al ( 2004). The three haloes 
studied, labelled A, B, C, are the most massive haloes in the simu- 
lation excluding haloes whose real density map shows an ongoing 
major merger. 

In each of the analysed haloes, we find «(/) to be well de- 
scribed by a power-law, 

-2.50±0.05 



vif) « r 



(37) 




over 3 to 5 decades in /. It is typically valid between slightly above 
/vir and slightly below f2o%- Outside this range, v{f) gradually 



Figure 4. The volume distribution of pha.se-space density, v(f), for each 
of the nine haloes analysed in this paper (see Table l4. iT . The curves were 
shifted to coincide at / = /-2, where the local log slope of is —2, and 
were then shifted vertically by 4 decades relative to each other. A power- 
law line f (/) oc f~'^'^ is shown on top of each curve. Marked on each 
curve are the virial level /vir and the 20% statistical eiTor limit /2o%- The 
/ range is shorter for the L320 haloes because they were sampled with less 
particles. 



and systematically deviates downward from the power law. In the 
low-/ regime the deviation is associated with departure from the 
virial regime, and the high-/ deviation is consistent with being due 
to the limited mass resolution of the specific halo, as seen by the 
proximity to f2o% and as demonstrated in Appendix]^ The high-/ 
deviation from the power-law tends to occur at a smaller / value in 
L25, and even smaller in L320, due to the fact that Mvir is smaller 
respectively. 

There is no evidence for a significant dependence of v{f) on 
the halo mass. There may be a marginal trend for slight steepening 
of v{f ) as a function of mass, but only from steeper than /~^ '*^ at 
~ 10® M0 to flatter than /"^ '^^ at ~ 10^^ Mq. This indicates rel- 
ative insensitivity to the exact slope of the initial fluctuation power 
spectrum, which varies across the range from dwarf galaxies to 
clusters of galaxies. Additionally, the fact that we obtained essen- 
tially the same v{f) from simulations using two different numer- 
ical codes, indicates that the shape of v{f) is not an artifact of a 
particular simulation technique. 
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Table 1. Global properties of the 9 haloes analysed in this paper. The v(f) of each halo was calculated from all the A^cut particles that lie within a radius -Rcut- 
-Rvir and A/vir ths th2 virial radius and mass. C is the concentration of the halo, calculated from an NFW fit. /2 is the / level where the logarithmic slope of 
'v{f) is —2. This value, together with v{f2), are used to scale the v(f) curves of the different haloes in Fig.|4] Finally, "Code", z, rripar, '"force and erg describe 
the computer code, red-shift, mass-resolution, force-resolution, and normalisation of each simulation. 



Halo 


A^cut 




i?vir 




C 


f2^ 




Code 


z 


'^par 


^force 


0-8 






[kpc] 


[kpc] 


[Mq] 












[Mq] 


[kpc] 




LIa 


1.6 X 10*^ 


23 


21 


1.1 X 10"' 


7.0 


4.7 X 10* 


1.0 X 10-* 


ART 


2.33 


7.0 X 10^ 


8.7 X 10-2 


0.75 


Lis 


1.3 X 10® 


23 


20 


8.6 X 10^ 


4.3 


3.4 X 10** 


1.3 X 10"* 


ART 


2.33 


7.0 X 103 


8.7 X 10-2 


0.75 


Lie 


1.1 X 10® 


20 


19 


7.8 X 10^ 


7.5 


3.5 X 10* 


1.1 X 10"* 


ART 


2.33 


7.0 X 10^ 


8.7 X 10-2 


0.75 


L2bB 


1.1 X 10'' 


420 


320 


1.9 X 10^2 


17.4 


5.8 X 10^ 


1.3 X 10° 


ART 





1.2 X 10'^ 


1.4 X 10-1 


0.9 


L25c 


1.2 X 10" 


420 


330 


2.0 X 10^2 


12.8 


3.4 X 10^ 


4.0 X 10'' 


ART 





1.2 X 10'^ 


1.4 X 10-1 


0.9 


L25d 


1.4 X 10" 


420 


340 


2.2 X 10^2 


11.7 


3.4 X 10^ 


4.3 X 10'' 


ART 





1.2 X 10^ 


1.4 X 10-1 


0.9 


L320a 


4.6 X 10^ 


3,000 


2,700 


1.1 X 10^5 


6.28 


8.8 X 10^ 


3.4 X 10* 


TPM 





2.6 X 10'' 


4.7 X 10° 


0.95 


L320s 


4.5 X 10^ 


3,000 


2,700 


1.1 X 10^5 


5.00 


6.0 X 10^ 


6.7 X 10* 


TPM 





2.6 X 10^ 


4.7 X 10° 


0.95 


L320c 


4.6 X 10^ 


3,000 


2,700 


1.1 X 10^5 


6.43 


1.0 X 10^ 


2.4 X 10* 


TPM 





2.6 X 10« 


4.7 X 10° 


0.95 



units: A/qMpc ^km ^ 
units: Af-iMpc^'km^s-^ 



5 SUBSTRUCTURES 

5.1 dumpiness in Phase-Space Density 

Had /(x, v) been a function of tlie energy alone, and the haloes 
were completely spherical and isotropic, the power-law t;(/) oc 
/-2'^ would have implied via Eq. <25> that the real-space den- 
sity profile must also be a power law, in fact an isothermal sphere 



We see that the real-space density maps are dominated by 
the familiar relatively smooth trend of density decreasing from 
the centre outward, with several tight clumps spread through- 
out the halo. The phase-space density maps, on the other hand, 
are qualitatively different. While the global trend with radius is 
much less apparent, the subhaloes become the highest peaks, 
especially in the outer regions of the halo. For example, the 



p(r) oc r-2, at least over some finite range in r. This is clearly clumps with / > 10 A/q Mpc (fcm/s) (yellow-reddish 



not the case, as the simulated haloes are well described by a univer- 
sal average density profile whose local logarithmic slope is varying 
continuously from —3 at the outer parts to —1 or even flatter in the 
inner parts (SQ. We conclude that / is far from being a function of 
energy alone, and in particular the system must deviate significantly 
from spherical symmetry or isotropy. This could be mostly due to 
the clumpy substructure of the halo, where the surviving subhaloes 
contribute volume of high phase-space density to thus mak- 
ing it shallower than expected from a smooth system with an inner 
density slope flatter than —2. 

In order to address the hypothesis that is driven by sub- 
structure, we plot real-density and phase-space density maps of 
each halo in real-space slices. Fig.|5|and Fig. |6| show such maps 
for dwarf haloes B and C (all other haloes show a similar qualita- 
tive behaviour). 

The real-space density p of each particle was calculate d usin g 
a three-dimensional Voronoi tessellation ( van de Wevgaer t]ll994) . 
generated using the free qhull software package. We chose this 
technique to estimate the real-space density because it is very simi- 
lar in its adaptive nature to the Delaunay tessellation technique used 
to estimate the six-dimensional phase-space density. A brief ac- 
count of the Voronoi tessellation technique is found in AppendixIXI 
where it is used to estimate the errors in the DTFE method. 

The maps were produced in the following way. For each halo, 
we consider all the particles within an equatorial slice parallel to 
the ly-plane, whose width is 40% of the virial radius. The slice is 
divided into 500 x 500 x equal cubic cells, with n-^ set to have 
the cells cubic. The density (p or /) assigned to each cell is the 
average of the densities of all particle within it. ^From each group 
of cells with the same x and y, we plot the one with the highest 
density. 



colours) are found everywhere. The very high peaks, with / > 
lO^Af© Mpc-^(fcm/s)-^ (bright yellow colours), are all found 
at a considerable distances from the centre. The central peak in / 
is quite modest in comparison; the elongated structures near the 
centre of dwarf-halo B are most likely subhaloes in the process of 
merging. 

Fig. |7| highlights the same effect by showing p and / asso- 
ciated with a random subset of the N-body particles as a function 
of their distance r from the halo centre. A large portion of these 
particles follow the global trend of decreasing density with radius 
- they could be associated with a smooth-background component, 
for which / is approximately a function of energy alone. At radii 
r > 1 kpc, the high-/ values come in "spikes" corresponding to 
the subhaloes. While the spikes in p reach values comparable to 
the central peak, the spikes in / could be more than 100 times 
higher, indicating that the subhaloes are both compact and cold. 
We note in dwarf halo B, for example, that all the points where 
/ > IOI^A/q Mpc-^km-^s^ are in subhaloes. 

The other interesting feature of the spikes is the fact that they 
seem to get lower and broader as they get closer to the halo centre, 
and that beyond a certain radius (of about 2 kpc), the spikes com- 
pletely blend into the smooth background. This indicates that the 
subhaloes phase-mix and lose their high phase-space densities as 
they approach the halo centre. This seems to be the natural result 
of mergers and tidal effects, which both puff up the subhaloes and 
heat them up. 



5.2 Toy Model: Adding Up Small Haloes 

As a first attempt at trying to understand how the power law v{f) oc 
/-2'^ in a halo of mass M may originate from its substructure, we 
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45kpc 



Figure 5. Density maps of dwarf halo Llg in a slice of tliickness 0.4Rvir- 
Top: real-space density. Bottom: phase-space density. The units in the 
colour key are log(p/[A/Q Afpc^'^]) and log(//[AfQ Afpc^''km~''s^]) 
respectively. The very-high / values are found inside clumps which are 
typically far away from the halo centre. 



simply add up the typical contributions from the general popula- 
tion of haloes of different masses m smaller than A/, as predicted 
in the ACDM cosmolog y. Based on cos mological N-body sim- 
ulatio ns tiMoore et al . .19993 : Ichigna el al .2000. : .De Lucia et al. I 



|2004|), and in accordance with the Press-Schechter approximation 
lpress & Schechter ' '1974 and its extensions (e.g., iLacev & Cole I 
[1993; Sheth 2003), the mass function of small-mass haloes (not 
necessarily subhaloes) can be approximated by dn/dm oc m""^, 
with 7 ~ 1.8 — 2.0. Additionally, we assume that the average den- 
sity profiles of haloes of different masses have the same functional 
form and are simply scaled versions of each other. This is estab- 
lished by n-body simulations for the case of i solated haloes, but i s 
less clear when one considers subhaloes (e.g., lHavashietatl2003h . 
Nevertheless, we adopt this assumption as our first crude toy model. 

As a first approximation we assume that the haloes all form at 
the same time in an Einstein-deSitter cosmology, so they have the 




45 kpc 



Figure 6. Density maps of dwarf halo Lie- See Fig.lsl 



same characteristic real-space density pm — p (a. constant factor 
times the universal density), and therefore their typical radii scale 
like Tm oc m^^"^. Based on the virial theorem, the velocity disper- 
sions then scale like cr„i oc m^^'^. Therefore, the typical phase- 
space volume of a halo of mass m scales like Vm oc r^Cm oc m^, 
and its typical phase-space density is fm oc m/Vm oc m~^. If we 
denote by v{f) the universal, un-scaled, normalised, dimensionless 
probability distribution function relevant for all the haloes, then the 
volume distribution function Vm{J) of a halo of mass m, defined 
such that Vm{f)df is a volume, is given by 



Jm V Jm ' 



(38) 



Then the total contribution to of halo M from the population 
of smaller haloes m < \iM is 

'"if) ^ / -;—Vm{f)dm^ / m^^^v{mf)dm, (39) 
Jo ^"^ Jo 
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oc m(4+-'/2, oc m-(2+-)/2, and we finally 



Figure 7. Densities as a function of radius for dwarf halo Lis, using a ran- 
dom set of 4% of the particles. Top: real-space density. Bottom: phase-space 
density. The background particles define a general trend of decreasing den- 
sity with radius, while the spikes con'espond to subhaloes. The phase-space 
density spikes are higher than the central peak because the subhaloes are 
cold. They become shorter and broader at smaller radii, indicating heating 
and puffing up by tidal effects and mergers. 



where [i,M is the mass of the largest subhalo (/i < 0.5). Changing 
variables m mf, we finally get 



vif) = / 



-(4-7) 



^v{s) ds . 



(40) 



Thus, v{f) is a multiple of a power-law and a monotoni- 

cally increasing function of /; it therefore has to be shallower than 
the power law For 7 ^ 1.8, this means that v(f) is shal- 

lower than f~^'^, which is significantly shallower than the mea- 



sured / 



. This idealised toy model does not seem to work. 



This calculation can be generalised to the case where haloes 
of smaller masses form first, as implied from the slope of the initial 
fluctuation power spectrum Pk- For an Einstein-deSitter cosmol- 
ogy, the formation time scales with m such that pm oc m^" , with 
V = (n + 3)/2, and n is the local power index, n — din Pk/dln k, 
at the relevant effective scale for masses < /xAf. Then, in anal- 
ogy to the calculation in the previous paragraph, we have am oc 



m 

obtain that v{f) must be flatter than the power law 

V(f) oc /-2-2(2-7)/(2+.) 



(41) 



For n — —3, namely 1/ = 0, the asymptotic index for small haloes 
where all haloes form at the same epoch, we recover the earlier 
result that v{f) is flatter than For ACDM on galactic 

scales, the effective power index is n ^ —2.3, so ^ 0.35, and 
for 7 ^ 1.8 we obtain that v{f) is flatter than /~^'^^. We see that 
the time dependence makes only a little difference. 

The total v{f) of the halo should be the sum of the background 
contribution and the subhalo contribution, but at the high-/ range 
we expect v{f) to be dominated by the contribution of subhaloes, 
as seen earlier in this section. The above toy model thus predicts 
that should be flatter than /~^'^, in conflict with the measured 

vif) cc r^'. 

We conclude that a halo is not simply an ensemble of clumps 
drawn from the general population of smaller haloes. The sub- 
haloes may have a different mass function, their shape properties 
may vary differently with mass, and they both could vary with dis- 
tance r from the host-halo centre. If one keeps the scaling relation 
with = and ignores any variation with r, the required subhalo 
mass function for matching the measured v{f) oc ,/ ^^'^ is with 
7 = 1.5 (compared to 7 = 1.8). Indeed, a flattening of the subhalo 
mass function could be a natural result of the inevitable dynamical 
evolution of the subhaloes in the potential well of their host halo. 
The phase mixing due to tidal effects, including total disruption, 
is likely to be more efficient in less massive subhaloes, thus flat- 
tening the mass function. Also, dynamical friction is more efficient 
in making more massive subhaloes sink into the halo, thus making 
the mass function flatter in the inner parts. However, recent sim- 
ulations i ndicate that the subha lo mass function is not flatter than 
7 ~ 1.7 jPe Lucia et al. |2004) . indicating that the tidal effects on 
the inner structure of subhaloes must also have an important role. 
These are matters for more detailed future studies, but the failure 
of the idealised toy model analysed above to reproduce the magic 
power law v{f) oc f~^'^ indicates that the phase-space density is 
likely to provide a useful tool for studying the dynamical evolution 
of subhaloes in host haloes. 



6 DISCUSSION AND CONCLUSION 

Using Delaunay tessellation, we developed a method for measur- 
ing the 6-dimensional coarse-grained phase-space density /(x, v) 
in A'^-body systems. We focused, in particular, on measuring the 
phase-space volume distribution function, v{f). We applied this 
technique to several simulated haloes of ~ 10^ particles, formed 
by hierarchical clustering in the standard ACDM scenario, and ob- 
tained two striking new results. 

First, v{f) is well described by a power law, v{f) oc 
j-2.5±o.05^ over 3 to 5 decades in /. The power-law regime starts 
at an / value which corresponds to the characteristic size of the 
virialized halo. It ends at an / value which is determined by the 
dynamical resolution limit of the specific simulation. Therefore, 
the true power-law range may extend to / ^ 00. This power law 
seems to be insensitive to the halo mass in the range 10^ — IO^^A/q, 
indicating insensitivity to the exact slope of the fluctuation power 
spectrum, as long as the haloes are built by hierarchical merging of 
clumps bottom up. 

Second, this power-law originates from substructures within 
the halo rather then the overall trend with radius. The substructure 
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completely dominates the high-/ parts of the v{f) distribution. The 
infalling clumps seem to phase-mix — by puffing up, heating and 
stripping — as their orbits decay from the virial radius inwards 
toward the halo centre and they melt into the halo smooth back- 
ground. 

Our first worry is that these results could be numerical arti- 
facts, or severely contaminated by such. Based on our error analysis 
and tests with mock datasets, we argue that the v[f) measured by 
the DTFE algorithm genuinely reflects the true phase-space prop- 
erties of the given A'^-body system over a broad range of /. The 
question is whether the phase mixing suffered by the subclumps as 
they approach the halo centre might be an artifact of numerical ef- 
fects such as two-body relaxation, leading to underestimated inner 
densities and/or overestimated internal velocities. A simil ar effect 
has b een pointed out using a one-dimensional toy model iBinnev I 
I2OO3I) . The apparent agreement between simulations run with dif- 
ferent codes and different resolutions is encouraging. In order to 
specifically address the effect of few-body relaxation, we intend to 
run twice a simulation of the same halo with the same number of 
particles but with a different force resolution (ongoing work with 
F. Stoehr). 

Assuming that the simulations genuinely reflect the true phys- 
ical behaviour under the Vlasov equation, the origin of the robust 
power-law shape of v[f) from the merging substructure becomes a 
very interesting theoretical issue. As demonstrated in ij5] a simple 
model using the mass function and the scaled profiles of the gen- 
eral halo population in the ACDM scenario does not reproduce the 
correct power law. This, and the apparent trend of the / spikes with 
radius, indicate that the structural and kinematic evolution of the 
subhaloes in the parent halo are important. Studies of tidal heating 
and stripping may be found useful in this modelling. 

It would be interesting to follow the phase-space evolution and 
the contribution to the overall v{f) by a single, highly resolved 
subhalo, or many of those, as they orbit within the parent halo and 
approach its centre. This may help us understand the nature of the 
interaction between the parent halo and its subhaloes, and the origin 
of the v[f) power law (ongoing works with E. Hayashi and with B. 
Moore). 

Another more general but speculative possibility is that the 
/"^■^ power law represents some sort of a cascade of relaxation 
processed in phase-space, in which high phase-space densities turn 
into lower (coarse-grained) densities through the process of mix- 
ing. In general, the fact that our findings are expressed in terms of 
the fundamental concept of phase-space density should make them 
more directly accessible to analytical treatment. In this respect, it 
may prove beneficial to investigate more closely the time evolu- 
tion of the v{f) of a cosmological halo and its components. This 
may shed light on the connection between the «(/) power-law be- 
haviour and the relaxation processes within the halo. 

We saw that the power-law behavior of v{f) is limited to 
the virial regime. It would be interesting to learn how this shape 
evolves in time as the halo virializes. A preliminary study (to 
be concluded and reported in another paper) indicates that in the 
intermediate-/ regime the «(/) of a pre-virialized system is sig- 
nificantly flatter than /~^'^, while in the high-/ regime it drops in 
a much steeper way. The }~^'^ behaviour seems to be a feature 
unique to virialized systems. 

We learnt that in the haloes that are built by hierarchical clus- 
tering, the power-law behaviour v{f) oc /^^'^ reflects the halo 
substructure. It would be interesting to find out whether this power- 
law behaviour actually requires substructure, or it is a more gen- 
eral phenomenon of virialized gravitating systems, valid indepen- 



dently of substructure. One way to answer this question would be to 
analyse simulated haloes in which all fluctuations of wavelengths 
smaller than the halo scale were removed, resulting in a smooth 
halo formed by monolithic collapse, with no apparent substruc- 
ture in the final configuration. As described in such haloes are 
known to still have NFW-like density profiles in real space, and one 
wonders whether they also have the magic power-law v{f). There 
are preliminary indications for a steeper v[f) in this case (Arad, 
Dekel & Moore, in preparation). If confirmed, it would indicate 
that the /^^'^ behaviour, while insensitive to the exact slope of the 
initial power spectrum, is unique to the hierarchical clustering pro- 
cess, and is not a general result of violent relaxation. 

Our current results are just first hints from what seems to be 
a promising rich new tool for analysing the dynamics and structure 
of virialized gravitating systems. The analysis could become even 
more interesting when applied to haloes including the associated 
gaseous and stellar components. 
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APPENDIX A: ERRORS IN Udei(/) 

The "true", underlying coarse-grained /(x, v) is what one would 
have measured using a fixed smoothing window in phase space (e.g. 
counts in fixed cells) and an arbitrarily large number of particles. 
Instead, the Delaunay Tessellation Rield Estimator uses adaptive 
cells in order to deal with the mass resolution limitations. There- 
fore, the relation between the measured Udoi(,f) and the underlying 
''truc(/) is not trivial. Both statistical and systematic errors might 
influence our results. We estimate these errors empirically below, 
but we start with an approximate model that provides a simple un- 
derstanding of the origin of the errors. 



Al A Voronoi Model: Fixed Smoothing v{f) 

A technique similar to the Delaunay tessellation, but somewhat 
simpler to interpret, is the Voronoi tessellation I van de Wevaae^ 
1994). Ror each particle, the Voronoi cell is defined as the region 
of phase space in which every point is closer to that particle than 
to any other particle. In this case, N particles define exactly A'^ 
Voronoi cells which cover all of phase space with no overlaps. If 
Vi is the Voronoi cell of particle i, then a natural mass-preserving 
way of estimating the phase-space density inside that cell is by 
fi = mi/\Vi\. We denote the quantities measured this way by /vor 
and Wvor (,/■). 

Much like the Delaunay tessellation, the Voronoi tessellation 
is an adaptive grid that enables one to estimate /(x, v) even in the 
presence of a relatively small number of particles. We use the De- 
launay method in our main analysis because it is somewhat more 
accurate (Schaao & van de Wevaaert 2000), and is easier to calcu- 
late numerically. However, the similar Voronoi method provides a 
simple way of learning about the properties of the measured v{f) 
and understanding the uncertainties associated with such a mea- 
surement. The empirical tests of the Delaunay measurements be- 
low demonstrate the relevance of the wisdom gained by analysing 
the Voronoi model. 

To understand the errors in the Voronoi density estimation, let 
us start with the trivial case where all of infinite phase-space is uni- 
formly filled with phase-space density /o, which is represented by 
an infinite number of particles with mass m. A volume V of phase- 
space would then contain on average a finite number of particles, 
foV/m. The Voronoi estimate of fi for each particle would fluc- 
tuate about /o due to the discreteness of the particle distribution. 
Since there is no typical scale in the problem, and each cell always 
contains one particle, the fluctuations Sf / f per Voronoi cell would 
remain at the same level even if one increases the average number 
density of particles while decreasing the mass of each particle in 
proportion, keeping /o the same. Therefore, the probability that the 
Voronoi estimated / would lie in the interval f f + df may 
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be written in terms of a universal probability distribution function 

Poo(///o)d(///o). 

When we consider a finite system in a box of volume V, with 
a finite number of particles A'^ inside it, we may expect the proba- 
bility distribution pjv ( ■ ) to deviate from its asymptotic form Poo ( • ) • 
As decreases, we expect pjv (•) to widen due to the increasing ef- 
fect of the boundaries. Next, examine a system with a non-uniform 
phase-space density, such as a cosmological halo. If the number of 
particles that represent this / is sufficiently large, we may approxi- 
mate every region in phase-space as being locally uniform, and es- 
timate its Voronoi / using the asymptotic Poc(-)- This assumption 
is expected to break down in regions with very high phase-space 
density, where the sampling may become poor and insufficient, or 
in regions where / has large gradients over small scales. Neverthe- 
less, lets assume for the moment that there exists an effective p(-) 
[not necessarily poo(-)] which properly approximates the fluctua- 
tion distribution of the Voronoi / for all particles. 

This assumption allows us to calculate the expectation value 
of Vvorif) for a finite system with a given /(x, v). This is done by 
first calculating the average of Vvor{f), the cumulative version of 
Vvorif) defined in J9}, and then differentiate it to obtain the average 
of Vvor(./). Assuming that /(x, v) is realised by 3> 1 particles, 
we divide phase space into a large number of cells uia, which are 
small enough to guarantee that (a) each cell is very unlikely to con- 
tain more than one of the TV particles, and (b) the value of /(x, v) 
is approximately constant within every cell. 

For each cell uia, we calculate the average of the 

contribution of this cell to Vvor(/). The contribution Va{f) would 
be non-zero only when there is a particle in the cell uia, a parti- 
cle whose assigned Voronoi estimated density is /vor- If fa is the 
true phase-space density in that the cell, then, according to our as- 
sumption, /vor would be chosen at random out of the probability 
distribution p(/vor//a) d.{fvoi/ fa) - Once /vor is chosen, Va{f) is 
given by 



Vaif) = 



m//v, 




(Al) 



Using our assumption (a) above, the probability of cja to host one 
particle is Pa = NM'^ /(x,v)dxdv, with M = Nm the 
total mass of the system. Therefore, 

/oo 
P(.fvor//a)d(/vor//a) m//vor (A2) 

/OO 
P[/vor//(x, v)]/;;id/vor , (A3) 

and so, 

(Kor(/)) = ^{14(/)) 

i 

' poo 

dxdv J p[/vor//(x,v)]/„"o^^d/vor 

; poo 

df'v(f') p(/vor//')/."ord/vor . (A4) 



In the last equality, we have used the exact v{f) to replace the six- 
dimensional phase-space integration. By differentiating Eq. <A4t 
with respect to /o, we finally obtain the desired expectation value: 



Kor(/))= / V{f')f'~'p{f/f')df' 

lo 



(A5) 



act v{f) with a fixed window function p{f /fft,^a)- The narrower 
p(-) is, the closer (wvor(/)) would be to the true «(/). However, 
as argued above, even when N —> oo the window p(-) does not 
approach a Dirac delta function; it rather converges to some finite- 
width distribution Poc.(-)- Therefore, even with an infinite resolu- 
tion the Voronoi tessellation would not produce the exact v{f)\ it 
converges to a convolution of it with a fixed window Poo(-)- 

This convolution would not affect the shape of the measured 
(vvorif)), and would preserve the true shape of v{f), provided 
that v{f) does not vary drastically over / scales which are smaller 
than the width of the window. In particular, when v{f) is a power 
law, the Voronoi algorithm would recover the same power law for 

{^^vor(/o))- 

A2 Empirical Testing with Mock Systems 

For an empirical study of the errors in the DTFE-measured v{f), 
and for testing how well the Voronoi model approximates these er- 
rors, we have performed a series of i'dci(/) measurements on sys- 
tems with known phase-space densities of the form 



-3/2 _„2/2<^2(^) 



(A6) 



corresponding to a spherical system in real space with a 
Maxwellian velocity dispersion. We have examined six such sys- 
tems with three different density profiles parametrised by a. 



Pa{x) 



a = 0, 0.5, 1.0 . 



and the following two types of dispersion profiles: 

/I// I T" ^ ^ 

Ov{x) 



M{x) 



crc(x) = 0.1 



(A7) 

(A8) 
(A9) 



The subscripts "v" or "c" denote a varying dispersion profile versus 
a constant one respectively. The «(/) for such systems is 



vif) = 
with 



:(/) 



X a {x)i 2 log 



dx . 



p{x) 



\2-Ka'^{x] 



n3/2 



(AlO) 



(All) 



We see that the measured (nvor (/)) is a convolution of the ex- 



and x{f) its inverse. 

Figure lATI shows the cumulative distributions of /del //true in 
different bins {fj,fj+i) of /true, for the a = 1 system and 
the Q = 1 cTc system. Both systems were realised using 10^ parti- 
cles. In both cases we have also plotted the cumulative distribution 
/del //true for a homogeneous Poisson distribution, realised in a 
six-dimensional cubic box with 10^ particles. We have verified that 
this distribution is essentially unchanged when the calculation is 
done with 10^ particles. Therefore, it should be regarded as good 
approximation to the asymptotic limit we would have reached in 
the different bins, had we used an infinite number of particles. The 
other four systems give essentially the same results. 

In the 1-v system, five bins were defined by fj — 
10''', 10"^ 10"^ 10"\ 10, 10^ We see that the shape of the dis- 
tribution, corresponding to the width of the differential distribution, 
is very much independent of /true- For all bins, the full width is 
less than a 1.5 decades. On the other hand, there does seem to be a 
systematic shift of the median toward larger values of /del //true 
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Figure Al. The cumulative distribution of /del //true as measured in dif- 
ferent bins of /true for two mock systems described by Eq. IA6I with 10® 
particles. Top: an NFW p{x) given by Eq. IA7I {a = 1.0) with a varying 
= M{x)/x. Bottom: same p{x) but with a constant crc{x) = 0.1. 
The solid black line represents the cumulative distribution of /del//truo for 
a homogeneous Poisson distribution, which is realised in a six-dimensional 
cube using 10® particles. This should be regarded as the asymptotic limit 
with an infinite resolution. 



for smaller values of /true- The large shift of the highest-/ bin 
(10 < / < 10"^) toward small /dei/,/true is a result of the very 
low number of particles in that bin, less than 500, which is insuffi- 
cient for representing such high phase-space densities. 

The bias toward larger phase-space densities in the low-J/j^.^^^ 
bins may be attributed to boundary effects: while the total mass 
of the N particles in the realisation is equal to the total mass one 
would obtain from the exact /(x, v) integrated over the infinite 
phase-space, the total phase-space volume used by the DTFE to es- 
timate /(x, v) is finite. It is the smallest possible convex polygon 
containing all TV particles. Therefore, we may expect an overesti- 
mate of /, which would be more pronounced near the boundaries. 
Nevertheless, as we shall see by comparing v^ciif) to Vtmcif), on 
scales of a few-decades, this bias is rather meaningless. 

The /dei/,/true distributions of the 1-c system are essentially 
the same as the ones in the 1-v system. Here 4 bins were defined 
by fi = 10"^10"^10"^10^10^ The overaH shape of the 
plots changes very little from bin to bin, and its width is about a 
decade and a half. Additionally there exists the bias towards larger 
/dei//true fatlos as /true gets Smaller. 
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Figure A2. The measured «dcl(/) versus the true v{f) for the three ctc 
mock systems, each performed with three different resolutions. 



To see how well the DTFE reconstructs the true i'(/), we have 
measured the i'dei(/) of the six mock systems using 10", 10^ and 
10® particles for each system. Fieures lA2l IA3I present the results 
of these numerical experiments. We see that with the highest res- 
olution, of 10® particles, the recovery is excellent over a range of 
almost 7 decades in /. Systematic deviations begin at the high-/ 
end and the low-/ end. At both ends, the deviations appear at about 
one to two decades inward to the highest and lowest values of /del 
in that realisation. When the resolution is decreased, the /-range 
where i'dei(/) closely matches ftrue(,/) narrows gradually. 

As argued above, the systematic deviations at the low-/ end 
are probably a result of the finite phase-space volume occupied by 
the particles. As the number of particles is increased, a larger por- 
tion of phase-space is covered, enabling the reconstruction of lower 
values of /del. 

The high-/ systematic deviations can be qualitatively under- 
stood using the Voronoi model and its convolution formula <A5t . 
Since the DTFE uses a finite number of particles to recover the / 
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Figure A3. Same as Fig. lA2l but for the CTv mock systems. 



field, there must be an upper cutoff, /i , for the spectrum of / val- 
ues that the DTFE can produce in principle. As far as the DTFE is 
concerned, the true v[f) has an effective cutoff at /i. If we plug 
this "true" truncated v{f) into the convolution formula <A5> . as- 
sume a distribution function p(-) with a width of about 1 decade, 
and set /i ~ 50, we recover the qualitative behaviour of the mea- 
sured i'dci(/) for the A'^ — 10® mock samples with CTv and the three 
values of a. This simple picture is, however, an oversimplification 
because the high-/ region is also effected by statistical fluctuations 
due to the low number of particles there. 

The convolution formula can also explain the overestimation 
by the Vdc\{f) in the high-/ regions of the mock systems with 
(Jc. In these systems, the transition between the low-/ power-law 
and the steep high-/ decline is rather sharp; it occurs over a scale 
comparable to the width of the window function p(-). As a result, 
the i'dci(/) measurement at a high-/ includes contributions from 
the higher v[f) values at lower /. This becomes very apparent in 
the case with q = and (Jc, in which / has an upper bound of 
/ ~ 63. While the true «(/) vanishes for all / values larger than 
this limit, the DTFE-measured v{f) vanishes only at an / value 



that is an order-of-magnitude larger, due to the I-decade width of 
the assumed probability distribution function p(-). 



A3 Statistical Errors in «dci ( /) • 

The Voronoi model can also be used to estimate the statistical er- 
rors in Uvor(/)- Strictly speaking, Vvoi{f) is an ill-defined random 
variable, as it is measured by differentiating Vvor(/), which is a 
super-position of Heaviside step functions, and as such «vor(/) is 
a sum of Dirac delta functions. Much like a white-noise process, 
its variance is infinite. In practise, however, we always compute 
Vvorif) by differentiating a smoothed version of Vvor(/) (using a 
spline, for example). Therefore, we may expect the statistical error 
in Vvorif) to be comparable to the statistical error in V^or(/)- The 
latter can be estimated in a way similar to how we estimated the 
average of Vvor (/) • 

To calculate ([AKor(/)]') = (VlXf)) " (Vvor(/))', we 
can use the definition of Va (/) to write 

Assuming that Va{f) is independent of Vs(/), the cross terms 
would cancel out from <^ [Al/„or(/)]^), leaving us with the upper 
limit 



([AKo.(/)]')<5](Fi(/)) 



(A13) 



Using arguments similar to those used for calculating 
{Va,(/)), one can show that 



m 
1 



dxdv / p[/vor//(x, v)]/„ord/vo(A14) 

/OO 
P[/vor//(x, v)]/iror d/viA15) 



j{VM)) 



Therefore, 



([AKo.(/)]')< j{Ko.(/)) , 
and the relative error A(/) is given by 



(A16) 



(A17) 



(A 18) 



Plugging m = M /N into the formula above, where M is the total 
mass of the system and A'^ the total number of particles, we recover 
the common large-numbers limit A(/) oc l/\/iV. 

To check how good this estimate is for the DTFE-measured 
v{f), we have measured the v{f) of 100 realisations of the a = 1 
(Jv mock system with A'^ = 10^ and A'^ = lO** particles, and 30 
realisations with A'' = 10^ particles. From these measurements we 
computed the true relative error in «(/) [with respect to the average 
of the DTFE-measured v{f), not with respect to the exact v{f)}, 
and compared it to the prediction of Eq. <A18> . Figure IA4| shows 
the comparison for the three resolutions. We see that Eq. <A18> 
performs well as an upper bound for the statistical errors, except for 
the low-/ region. In that region the cumulative V{f) approaches a 
constant as / 0, due to the finite phase-space volume of the halo. 
This introduces fluctuations to v{f) as a result of the numerical 
differentiation in Eq. <33> . 

However, it is interesting to notice that whenever the statistical 
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Figure A4. The average (top) and relative error (bottom) in v[f) as mea- 
sured in tlie mock system (a = 1, Cv) sampled with three different res- 
olutions. The relative errors are compared to the analytic prediction of 
Eq. fAlsl . 

errors in Vdc\{f ) become important, they are overwhelmed by the 
low-/ or high-/ systematic errors. In that respect, the statistical 
errors in tidci(/) are of no big relevance. 



